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HIGHLIGHTS 


►  We  perform  pore-scale  simulations  for  different  XCT-reconstructed  electrode  morphologies. 

►  We  predict  cell  performance  for  different  flow  rates. 

►  Cell  voltage  saturates  at  high  flow  rates  when  concentration  approaches  uniformity. 

►  Fuel  starvation  is  found  to  have  a  negative  impact  on  cell  voltage. 

►  We  report  statistical  distributions  of  local  parameters  across  the  3D  porous  electrode. 
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A  3D  pore-scale  transport  resolved  model  is  used  to  study  the  performance  characteristics  of  a  vanadium 
redox  flow  battery  (VRFB)  with  various  electrode  morphologies  under  different  operating  conditions. 
Three  electrode  structures  are  reconstructed  from  X-ray  computed  tomography  (XCT)  images  of  porous 
carbon  felt  electrode  materials.  The  local  vanadium  concentration,  overpotential,  current  density  and 
overall  cell  voltage  for  the  positive  half  cell  are  examined.  The  results  indicate  that  the  cell  voltage 
increases  with  increasing  electrolyte  flow  rate  due  to  decreasing  concentration  gradients  of  vanadium 
species  within  the  porous  electrode.  However,  the  marginal  gain  in  cell  voltage  diminishes  once  the 
concentration  field  approaches  uniformity  under  convection-dominated  mass  transport  conditions  at 
sufficiently  high  electrolyte  flow  rates.  The  model  also  predicts  that  electrode  structures  with  low 
porosity  (high  surface  area)  result  in  more  uniform  and  lower  absolute  current  density  and  overpotential 
fields  at  the  expense  of  increased  pressure  drop.  Finally,  poor  cell  performance  is  observed  for  simula¬ 
tions  operated  at  low  electrolyte  flow  rates  and  low  states  of  charge  due  to  the  fuel  starvation  (i.e., 
insufficient  amount  of  reactant  in  the  cell). 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  vanadium  redox  flow  battery  (VRFB)  is  considered  a  high- 
capacity  energy  storage  technology  used  for  storing  energy  from 
intermittent  energy  sources  such  as  wind  and  solar  [1—3].  At  its 
core,  the  VRFB  is  a  rechargeable  battery  that  consists  of  two  porous 
electrodes  composed  of  carbon  fibers  that  form  an  electrically 
conductive  fibrous  network  separated  by  an  ion  exchange 
membrane  (shown  in  Fig.  1).  Vanadium-based  electrolytes  are 
circulated  from  external  storage  tanks  through  the  porous  elec¬ 
trodes,  where  electrochemical  reactions  occur  on  the  surface  of  the 
electrode  fibers.  The  vanadium  species  are  dissolved  in  a  sulphuric 
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acid  (H2SO4)  solution  where  V3+  and  V2+  ions  are  present  in  the 
negative  electrolyte,  and  VOj  and  V02+  (also  known  as  V5+  and 
V4+,  respectively)  ions  are  present  in  the  positive  electrolyte. 
During  the  discharging  cycle,  the  following  reactions  take  place  at 
the  surface  of  the  carbon  fibers: 

Negative  half  cell  :  V2+  -►  V3+  +  e_  (1 ) 

Positive  half  cell  :  VOj  +  2H+  +  e“  ^  V02+  +  H20  (2) 

The  open  circuit  voltage  (OCV)  for  reaction  (1 )  is  E°~  =  -0.255  V 
[4]  and  that  for  reaction  (2)  is  E0,+=0.991  V  [5]  based  on  the  Gibbs 
free  energy,  leading  to  a  theoretical  standard  cell  OCV  of  1.246  V  at 
a  temperature  T  =  298 K.  Here,  the  superscripts  -  and  +  denote  the 
negative  and  positive  half-cell  reactions,  respectively. 
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Fig.  1.  Schematic  of  the  positive  half  cell  of  a  vanadium  redox  flow  battery. 


The  main  advantage  of  VRFBs  is  that  power  generation  and 
energy  storage  are  decoupled,  such  that  the  energy  storage  capacity 
is  determined  by  the  size  of  the  electrolyte  tanks,  whereas  the 
power  rating  is  dictated  by  the  size  and  the  number  of  redox  cells 
[6].  However,  transport  losses  in  the  electrode  and  low  electrolyte 
utilization,  among  other  losses,  significantly  affect  the  performance 
of  VRFBs.  During  VRFB  operations,  both  electrode  microstructures 
(e.g.,  active  surface  area  and  porosity)  and  electrolyte  flow  condi¬ 
tions  play  important  roles  in  determining  the  effectiveness  of 
electrolyte  transport  inside  a  VRFB  [7].  Experimental  identification 
of  optimal  electrode  microstructures  and  electrolyte  flow  configu¬ 
rations  for  VRFBs  is  very  challenging  and  expensive.  The  wide¬ 
spread  availability  of  supercomputers  has  made  the  use  of 
advanced  numerical  modeling  a  valuable  tool  in  assisting  the 
design  of  VRFBs. 

To  serve  this  purpose,  the  authors  have  recently  developed 
a  pore-scale  resolved  transport  model  [8]  for  the  VRFB  where  the 
porous  electrode  and  electrolyte  phases  are  distinctly  separated 
and  the  electrolyte  flow,  species/charge  transport,  and  electro¬ 
chemistry  are  coupled  and  solved  with  pore-scale  resolution.  The 
greatest  benefit  of  this  pore-scale  model  is  its  capability  of 
capturing  the  exact  electrode  morphology  and  electrolyte  flow 
conditions  and  how  these  parameters  affect  the  overall  cell 
performance.  As  opposed  to  volumetric  models  [9-15]  where  the 
electrode/electrolyte  domain  is  treated  as  a  single  continuum,  the 
pore-scale  resolved  model  does  not  rely  on  approximations  of  the 
electrode  microstructure  and  hence  offers  a  more  accurate  and 
detailed  understanding  of  the  coupled  electrolyte  flow,  species/ 
charge  transport  phenomena,  and  electrochemical  reactions 
occurring  locally  within  a  VRFB.  As  such,  the  pore-scale  model 


enables  investigation  of  key  structurally  dependent  parameters 
that  affect  the  electrolyte  flow  and  utilization,  thereby  allowing  for 
the  characterization  of  the  cell  voltage  parametrically  for  different 
electrode  structures  (e.g.,  size  and  orientation  of  the  fibers,  size/ 
distribution  of  pores)  and  different  operating  conditions  (e.g.,  flow 
rates,  electrolyte  concentration). 

In  this  paper,  the  pore-scale  model  developed  by  the  authors  [8] 
is  used  to  investigate  the  effects  of  electrode  structural  properties 
and  electrolyte  flow  conditions  on  the  electrolyte  utilization  and 
performance  characteristics  of  VRFBs.  A  3D,  isothermal,  positive 
half  cell  is  used  in  the  simulation,  where  the  electrode  geometry  is 
based  on  an  X-ray  computed  tomography  (XCT)  reconstructed 
porous  carbon  felt  electrode  material.  The  local  and  averaged 
vanadium  concentration,  overpotential,  current  density  and  cell 
voltage  are  reported  for  different  electrode  morphologies  (e.g., 
active  surface  area  and  porosity),  inlet  state-of-charge  (SOC),  and 
flow  rates. 

The  remainder  of  this  paper  is  organized  as  follows.  Section  2 
describes  the  XCT  reconstructed  electrode  geometry,  model 
assumptions  and  equations,  and  the  numerical  implementation  of 
the  pore-scale  resolved  model.  The  verification  of  the  half-cell 
model  is  shown  in  Section  3.  Results  are  then  presented  and  dis¬ 
cussed  in  Section  4.  Section  5  summarizes  conclusions  and 
recommendations  for  further  study. 

2.  Methodology 

2  A.  XCT-re constructed  electrode  geometry 

For  this  study,  a  sample  of  Electrolytica  GFS6-3  mm  carbon  felt 
(typically  used  in  VRFBs)  is  imaged  using  a  SkyScan  1172  X-ray 
tomograph.  A  virtual  volume  of  the  porous  electrode  medium  is 
then  reconstructed  from  the  segmented  tomograms.  The  3D 
geometry  is  finally  analyzed  for  key  structural  metrics  such  as 
porosity,  specific  surface  area,  pore  and  fiber  sizes,  and  phase 
connectivity.  A  thorough  discussion  of  this  methodology  is 
provided  in  Qiu  et  al.  [8]. 

Fig.  2  shows  three  different  electrode  structures  that  were 
reconstructed  from  XCT.  The  geometric  parameters  corresponding 
to  these  structures  are  given  in  Table  1.  The  increased  fiber  density 
structures  are  attained  by  compressing  the  physical  electrode 
medium,  a  phenomenon  that  occurs  during  the  assembly  of  actual 
VRFB  cells.  This  increase  in  fiber  density  also  results  in  lower 
porosities  and  higher  specific  surface  area.  The  structure  of  Fig.  2(a) 
is  from  an  uncompressed  sample  of  electrode  and  optimally 
matched  the  experimentally  determined  geometric  parameters  of 
the  virgin  electrode  medium.  All  of  the  electrode  structures  occupy 
a  volume  of  dimensions  900  x  135  x  450  pm3  (200  x  30  x  100 
voxels3  at  4.5  pm  per  voxel  resolution).  A  convention  of  dimen¬ 
sionless  length,  width,  and  height,  X*,  Y*,  Z*,  respectively,  are  used 
throughout  this  paper,  where  each  dimension  is  normalized  rela¬ 
tive  to  their  respective  total  distance,  i.e.  X*  =  X/L,  Y*  =  Y/Wf  and 
Z*  =  Z/H. 

2.2.  Model  equations 

Pore-level  transport  models  are  developed  to  simulate  the  flow 
of  electrolyte,  the  transport  of  ionic  species  (reactants  and  prod¬ 
ucts),  and  transport  of  charge  (electric  current)  through  both  the 
solid  and  liquid  phases,  as  well  as  the  electrochemical  reactions  at 
the  electrode  fiber  surface.  The  governing  equations  for  these 
phenomena  are  summarized  in  Table  2  and  explained  in  detail  in 
Qiu  et  al.  [8].  The  electrolyte  is  assumed  to  be  Newtonian  with 
a  uniform  viscosity  and  the  flow  is  incompressible,  driven  by 
a  constant  pressure  gradient  in  the  Z-direction,  and  with  a  no  slip 
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Fig.  2.  3D  model  geometries  that  are  used  in  the  pore-scale  simulation.  The  fiber  structures  are  digitally  reconstructed  from  XCT  imaging. 


boundary  condition  at  the  fiber  surface.  The  electrolyte  is  assumed 
to  be  dilute  for  all  ionic  species  and  the  sulfuric  acid  completely 
dissociates  into  sulfate  (SO4-)  and  protons  (H+).  The  model 
parameters  used  in  this  study  are  given  in  Table  3. 

In  the  present  study,  only  the  positive  half  cell  is  simulated.  As 
such,  only  the  transport  of  species  je  {V4+,  V5+,  H+,  H20}  is 
considered.  It  is  convenient  to  define  the  inlet  concentration  of  any 
species  (Cjn)  in  terms  of  the  total  vanadium  concentration  (C°,+) 
and  initial  proton  concentration  (C°v+)  with  respect  to  the  state-of- 
charge  (SOC)  of  the  electrolyte  using  the  following  expressions: 


Table  1 

Geometrical  properties  of  XCT-reconstructed  carbon-felt  electrode. 


Structure 

a. 

b. 

c. 

Porosity 

- 

0.932  (0.92603) 

0.897 

0.845 

Specific  surface  area 

m-1 

29,600  (37,500a) 

47,700 

69,700 

Mean  pore  diameter 

pm 

197  (102.2  ±  8.27a) 

122 

92.7 

Mean  fiber  diameter 

pm 

17.1  (20.80  ±  6.53b) 

16.7 

17.6 

a  Obtained  from  mercury  intrusion  porosimetry. 
b  Obtained  from  SEM  imaging. 


SOC  =  Cv/C°’+ 

c^n  =  c0+soc 

C&  =  C°’+-(l  -  SOC) 
Cj,n;+  =  C°v+  +  C0  +  •  SOC 

M  M 


(3) 


where  the  total  vanadium  and  initial  proton  concentrations  are 
given  in  Table  3.  Unless  otherwise  stated,  the  SOC  for  all  simulations 
is  50%. 

A  simplified  treatment  of  the  migration  of  H+  and  H2O  across 
the  ion-exchange  membrane  is  adopted  for  this  work  where  the 
bulk  generation  and  depletion  of  H+  and  H20  in  the  half  cell  are 
accounted  for  by  treating  them  as  participants  in  the  electro¬ 
chemical  reactions  at  the  active  surface  with  a  stoichiometric  rate 
of  molar  flux  according  to  Eq.  (2)  [8].  The  ion  exchange 
membrane  in  this  work  is  used  as  a  grounded  reference  electrode, 
as  seen  in  Fig.  1.  It  should  also  be  noted  that  the  migration  term 
accounting  for  the  effect  of  voltage  gradients  has  a  negligible 
effect  in  redox  flow  batteries  [15]  and  therefore  is  ignored  in  this 
study. 
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Table  2 

Governing  equations  and  boundary  conditions  [8]. 


Conservation  equations 


Boundary  conditions 


Fluid,  electrolyte 

Vu  =  0, 

aii 

—  +  (U-V)U  = 

1  9 

=  — Vp  +  yVzu 

P 

Co 

qjN 

C  II 

o 

=  Pin,P|z‘=l 
au 

_  pout 

=  0 

my 

=o  9y|y=i 

"1  fiber  =  0 

Species,  electrolyte 

hCj 

=  D;V* 2Q  +  V- 

Cjlz*=0  = 

.  rin 

'  j 

at  J 

J  J 

L  RT  J 

aqi 

=  aqi 

_8Cj| 

®z  lz*=1 

[  9y  y*=0 

1 

if 

Charge,  electrode 


Charge,  electrolyte 


Electrochemistry  (positive  half  cell) 

Reaction  rate 

Overpotential 


V-(KsV0)  =  0 


NIV-n  =  -Ny-n  =  k(Cfv)“c(q)“a  [exp(^j^)  -expf-^- 
V  =  0s-0e-[E°+(RTlF)\n(CvlCw)]  1  V  J  V 


dXly*- 


=  Jex  t) 


/Cc— —  i  /Cc-— —  i  K s~ ^ i  /Cs— —  —  0 

dy|y*=0  9y|y*=l  ®Z\z*=0  ®Z\z*=t 


[(§Y,zJDici)V(i>+FY,ziDi'*ci\  =° 

0|**=O  -  0 

a^l  d(f)\  801  801 

K-p— —  ,  K-p— —  ,  K-p— —  ,  K-p— — 

9y|y*=0  9y|y*=1  dz|z*=o  dz|; 

2.3.  Numerical  implementation 

The  numerical  program  is  written  in  FORTRAN  90  using  the 
Message  Passing  Interface  (MPI)  for  parallel  processing.  Simula¬ 
tions  are  initiated  by  supplying  the  numerical  model  with  the 
domain  geometry,  physical  parameters,  inlet  velocity,  SOC,  and 
external  current  density  at  the  current  collectors.  The  lattice 
Boltzmann  Method  (LBM)  [17]  is  then  employed  to  solve  for  the 
flow  field  of  the  electrolyte  inside  the  pore  space  of  the  electrode 
medium.  A  fully  implicit  3D  finite-volume  method  (FVM)  [18]  is 
then  used  to  solve  the  coupled  charge  and  species  transport 
equations  in  the  domain.  A  typical  simulation  of  the  XCT  structures 
with  a  relative  tolerance  of  1  x  10-8  for  all  variables  converges 
within  20  h  using  a  112-core  cluster  running  on  the  NSF  XSEDE 
system.  Details  of  the  numerical  implementation  are  provided  in 
Ref.  [8]  and  a  schematic  of  the  multiphase  simulation  process  is 
shown  in  Fig.  3. 

The  fluid  velocity  and  pressure  as  calculated  from  LBM  are  in 
lattice  units  and  are  related  to  physical  units  through  appropriate 
dimensionless  quantities.  As  such,  the  physical  fluid  velocity  is 
obtained  by  matching  the  Reynolds  number  (Re)  in  both  lattice  and 
physical  units  as  follows: 


b^avg^avg\  _  /f^avgf'avg 

v  J  LBMunits  V  v 

where  the  characteristic  length  scale,  Lavg,  is  taken  to  be  the  average 
pore  size  based  on  the  XCT  data,  the  characteristic  velocity,  Davg,  is 
defined  as  the  average  velocity  at  the  cell  inlet  (Z*  =  0),  and  the 
actual  control  volume  size  is  the  length  scale  (Ax  =  1)  in  LBM  units. 

A  prescribed  pressure  difference  between  the  inlet  and  outlet  is 
specified  such  that  the  desired  flow  rate  is  attained.  Following 
a  similar  process,  the  applied  pressure  across  the  inlet  and  outlet  of 
the  flow  cell  can  be  converted  from  lattice  units  to  physical  units 
from  the  Euler  number  by  using  the  following  expression: 


physicalunits 


(4) 


\A P_\ 

P^avg/  LBMunits 


physicalunits 


(5) 


3.  Half-cell  model  verification 

The  half-cell  configuration  of  the  present  study  is  a  modification 

of  the  full-cell  model  that  was  previously  validated  against  pub¬ 

lished  experimental  and  simulated  results  [8].  Further  verification 


Table  3 

Operating  parameters  used  in  the  simulation. 


Description 

Symbol 

Value 

Units 

Electrolyte  density 

P 

1000 

kg  ITT3 

Electrolyte  viscosity 

V 

1.0  x  10"6 

m2  s"1 

Solid  electrode  conductivity 

Ks 

1000a 

S  m1 

Total  vanadium  concentration 

C°’+ 

2000b 

mol  m  3 

Initial  proton  concentration 

f-0,+ 

4000b 

mol  m  3 

Initial  water  concentration 

r° 

*-h2o 

4200a 

mol  m  3 

Anodic  transfer  coefficient 

aa 

0.5a 

- 

Cathodic  transfer  coefficient 

ac 

0.5a 

- 

Standard  reaction  rate  constant 

k 

6.8  x  10"7  [16] 

m  s"1 

Equilibrium  potential 

E° 

0.991  [5] 

V 

Operating  temperature 

T 

298 

K 

External  current  density 

Jext 

400 

A  m  2 

Species  Symbol  (C,) 

Charge  (Zj)  Diffusivity  [m2  s  1] 

V02+  Civ 

+1 

3.9  x 

0 

0 

w 

VOJ  Cv 

+2 

3.9  x 

0 

0 

UJ 

H+  C£ 

+1 

9.312 

x  10-9  [15] 

H2O  Ch20 

0 

2.3  x 

10"9  [13] 

sol-  csor 

-2 

2.2  x 

0 

0 

UJ 

a  Estimated. 

b  Based  on  vanadium  solubility  limit. 
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Fig.  3.  Schematic  of  the  pore-scale  multiphase  simulation  geometry  shown  in  a  2D 
plane,  (a)  Solid  fibers  (green)  and  pore  spaces  (grey)  are  assembled  in  a  Cartesian 
lattice,  with  the  interface  between  the  two  phases  (bold  black  line)  signifying  the 
active  surface  where  electrochemical  reactions  take  place,  (b)  4  x  4  subsection 
showing  the  electrochemical  reaction  occurring  at  an  active  surface  for  discharge  of 
the  positive  half-cell,  which  is  modeled  using  the  Butler— Volmer  equation.  The  fluid 
transport  in  the  pore  phase  is  solved  using  the  LBM,  while  the  species  and  charge 
transport  is  solved  using  the  FVM.  (For  interpretation  of  the  references  to  color  in  this 
figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 
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Table  4 

Operating  parameters  (adopted  from  You  et  al.  [15])  for  the  verification  study. 


Description 

Symbol 

Value 

Units 

Porosity 

E 

0.929 

- 

Specific  surface  area 

Ae 

16,200 

m1 

Half-cell  width 

L 

0.003 

m 

Average  velocity 

Wave 

0.014 

ms"1 

Initial  vanadium  concentration 

C° 

1500 

mol  m  3 

Equilibrium  potential 

E° 

1.004 

V 

External  current  density 

Jext 

-400 

Am"2 

of  the  3D  half-cell  model  is  performed  in  this  section  against  the 
available  half-cell  simulation  results  that  are  presented  for  the  2D 
volumetric  model  of  You  et  al.  [15].  Table  4  summarizes  the 
parameters  that  are  adopted  for  this  verification  study.  Note  that 
a  specially  designed  electrode  structure  is  used  [8]  in  order  to 
preserve  the  geometric  parameters  used  in  the  2D  volumetric 
model.  Fig.  4(a)  shows  the  averaged  solid  electrode  potential  as 
a  function  of  the  inlet  SOC,  while  Fig.  4(b)  shows  the  averaged 
surface  overpotential  along  the  X-direction  for  the  positive  half  cell 
under  galvanostatic  charge  of  400  A  nrr2,  which  have  average 
relative  errors  of  0.38%  (3.9  mV)  and  10.2%  (0.41  mV),  respectively. 
The  slight  discrepancies  observed  are  most  likely  the  contributions 
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Fig.  4.  Comparison  of  simulation  results  obtained  using  the  3D  pore-scale  model  and 
the  2D  volumetric  model  of  You  et  al.  [15]  for  the  positive  half  cell  under  galvanostatic 
charge  of  400  A  m-2.  (a)  Averaged  solid  electrode  potential  (average  error  =  0.38%)  and 
(b)  averaged  X-direction  surface  overpotential  (average  error  =  10.2%)  at  50%  SOC. 


from  the  localized  effects  of  fluid,  species  and  charge  transport  that 
can  only  be  captured  using  the  pore-scale  model. 


4.  Results  and  discussion 

4.1.  Effect  of  electrolyte  flow  rate 

Flow  rate  is  an  important  control  parameter  in  the  operation  of 
the  VRFB  since  low  flow  rates  induce  poor  electrolytic  circulation 
and  yields  stagnant  regions  (in  terms  of  electrochemical  activity)  in 
the  electrode  while  flow  rates  that  are  too  high  pose  leakage  risks 
and  large  pumping  power  requirements.  In  this  study,  flow  rates 
that  correspond  to  Re  number  of  order  1CT5  to  10  are  simulated  for 
the  93.2%  porosity  structure  shown  in  Fig.  2(a)  under  galvanostatic 
discharge  at  400  A  m-2,  with  a  specified  inlet  SOC  of  50%.  For 
comparison,  the  Re’s  that  are  used  in  simulations  [9,10,12,13,19]  and 
experiments  [20]  range  from  10-1-1  and  10-2-10,  respectively. 

The  effect  of  flow  rate  on  the  total  half-cell  voltage  and  the 
associated  pressure  drop  is  shown  in  Fig.  5.  The  results  indicate 
that  performance  improves  with  increasing  flow  rate.  For  low 
flow  rate  cases  where  Re  <  0.005,  it  is  anticipated  that  mass 
transfer  processes  become  diffusion-dominated,  which  result  in 
large  gradients  in  the  concentration  field.  It  also  appears  that  for 
a  flow  rate  of  Re  =  0.05,  the  gain  in  cell  performance  with 
increased  flow  rate  becomes  negligible  under  the  prescribed 
operating  conditions.  Shah  et  al.  and  Al-Fetlawi  et  al. 
[9,10,12,13,19]  also  observed  that  marginal  gains  in  output  volt¬ 
ages  decrease  with  increasing  flow  rate,  although  their  study  used 
a  more  limited  range  of  flow  rate.  Our  results  indicate  that  there 
is  a  critical  flow  rate  (approximately  Re  =  0.05  under  these 
conditions)  above  which  the  output  voltage  saturates.  As  such, 
operating  the  VRFB  at  higher  flow  rates  does  not  justify  the 
exponentially  increasing  pumping  power  required. 

Fig.  6  presents  the  local  distributions  of  V5+  concentration 
normalized  relative  to  its  specified  inlet  value  (herein  referred  to 
as  concentration  or  electrolytic  concentration)  and  surface 
current  density  for  flow  rates  corresponding  to  Re  =  0.05,  0.005, 
and  0.0005.  From  Fig.  6(1)  it  appears  that  the  circulation  of 
species  improves  with  increasing  flow  rates.  The  large  concen¬ 
tration  and  current  density  gradients  in  the  diffusion-dominated 
cases  can  clearly  be  seen  in  Fig.  6(al)  and  Fig.  6(a2),  respec¬ 
tively.  By  contrast,  more  uniform  distributions  of  concentration 
are  observed  from  the  convection-dominated  mass  transfer 
conditions  in  Fig.  6(cl),  which  indicate  better  circulation  of  the 
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Fig.  5.  Effect  of  flow  rate  on  half-cell  voltage  and  applied  pressure  drop.  Simulations 
are  performed  on  the  93.2%  porosity  half-cell  electrode  subject  to  galvanostatic 
discharge  of  400  A  m-2  with  a  specified  inlet  SOC  of  50%. 
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Fig.  6.  Distribution  of  local  transport  properties  shown  on  (1)  a  slice  plot  for  the  normalized  V5+  concentration  (C/Cin  vs+ )  about  the  Y-midplane  (note  the  different  scale  bars),  (2) 
a  surface  plot  of  the  current  density  (/)  on  the  3D  electrode  fibers,  and  (3)  a  scatter  plot  histogram  of  the  corresponding  J  and  C/Cin  vs+  at  every  applicable  control  volume  on  the 
electrode  structure  (using  32  x  32  bins).  Results  are  shown  for  different  flow  rates  corresponding  to  Re’s  of  (a)  0.0005,  (b)  0.005,  (c)  0.05.  Simulations  are  performed  on  the  93.2% 
porosity  half-cell  electrode  subject  to  galvanostatic  discharge  of  400  A  m-2  with  a  specified  inlet  SOC  of  50%. 


reactants.  In  Fig.  6(3),  the  local  concentration  and  current  density 
at  every  discrete  active  surface  on  the  electrode  structure  are 
collected  in  a  scatter  plot  histogram.  It  can  be  seen  that  a  larger 
scatter  indicates  more  diffusion-dominated  mass  transfer  condi¬ 
tions  as  a  result  of  larger  concentration  and  current  density 
gradients.  The  lower  diagonal  slope  seen  in  Fig.  6(a3),  Fig.  6(b3), 
and  to  a  lesser  degree  in  Fig.  6(c3)  suggests  that  the  highest  rates 
of  reactions  occur  in  regions  with  the  highest  electrolytic 
concentration. 


Fig.  7  shows  the  spatial  profiles  of  concentration,  overpotential 
and  current  density  in  the  X  and  Z-direction  as  well  as  their 
statistical  distributions  for  all  six  simulated  flow  rates.  Note  that  for 
flow  rates  of  Re  >  0.05,  the  spatial  and  statistical  distributions  for  all 
three  parameters  are  nearly  indistinguishable  from  each  other, 
which  indicates  that  (i)  convection-dominated  mass  transfer 
conditions  appears  to  exist  at  this  flow  rate  and  (ii)  increasing  the 
flow  rate  over  the  critical  value  ( -  Re  =  0.05  under  these  operating 
conditions)  will  yield  limited  changes  to  the  spatial  distributions  of 
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Fig.  7.  Effects  of  flow  rate  on  averaged  distributions  of  (a)  normalized  V5+  concentration,  (b)  overpotential,  and  (c)  current  density  along  (1 )  the  Z-direction  and  (2)  the  X-direction. 
Statistical  information  for  each  parameter  in  the  entire  domain  is  shown  in  (3),  where  the  dots  indicate  the  mean  value,  the  capped  bars  indicate  the  variation  to  one  standard 
deviation,  and  the  uncapped  lines  indicate  the  minimum  and  maximum  values.  Simulations  are  performed  on  the  93.2%  porosity  half-cell  electrode  subject  to  galvanostatic 
discharge  of  400  A  m-2  with  a  specified  inlet  SOC  of  50%. 


concentration,  overpotential  and  current  density.  For  the  diffusion- 
dominated  cases,  the  high  propensity  for  incoming  reactants  to 
oxidize  near  the  inlet  manifests  itself  in  the  steep  gradients  in  all 
three  parameters  in  the  Z-direction,  as  seen  in  Fig.  7(1).  For  flow 
rates  of  Re  >  0.05,  gradients  in  current  density  and  overpotential 
are  most  pronounced  in  the  direction  perpendicular  to  the  current 
collector,  as  seen  in  Fig.  7(b2)  and  Fig.  7(c2).  It  can  be  observed  from 
Fig.  7(a2)  that  most  reactions  occur  near  the  current  collectors  due 
to  an  increase  of  effective  conductivity  in  the  liquid  electrolyte 
relative  to  the  solid  fibers  as  a  result  of  the  high  porosity. 

The  statistical  distributions  shown  in  Fig.  7(3)  are  evaluated  for 
every  active  surface  in  the  entire  3D  simulation  geometry.  Fig.  7(a3) 
shows  that  the  average  concentration  decreases  and  that  increasing 
concentration  gradients  are  present  at  lower  flow  rates.  Fig.  7(b3) 


and  Fig.  7(c3)  suggests  that:  (i)  the  average  overpotential  and 
current  densities  remain  constant  for  all  flow  rates,  as  required  by 
Faraday’s  law;  (ii)  variations  in  these  parameters  increase  with 
decreasing  flow  rate.  These  observations  indicate  that  cell  voltage 
will  increase  with  increasing  flow  rates  and  abruptly  saturate  when 
convection-dominated  mass  transfer  conditions  are  established 
(Re  >  0.05  under  these  operating  conditions).  Past  this  point,  the 
concentration  field  approaches  uniformity  (to  the  value  specified  at 
the  inlet),  and  so  open  circuit  voltages  saturate  along  with  the 
overall  cell  voltages.  It  can  be  expected  that  with  a  larger  simulation 
domain,  concentration  gradients  will  be  larger,  and  so  cell  voltage 
will  saturate  at  higher  flow  rates.  It  is  important  to  note  that  for  the 
diffusion-dominated  cases,  nonphysical  positive  values  for  over¬ 
potential  and  current  densities  are  observed.  These  are  purely 
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a  mathematical  artifact  of  the  logarithmic  term  in  the  Nernst 
potential  (which  can  be  very  large  in  areas  of  low  concentration), 
and  are  found  to  occupy  less  than  1%  of  the  entire  total  active 
surface. 

4.2.  Effect  of  electrode  structure 

The  performance  of  VRFBs  depends  strongly  on  the 
morphology  of  the  electrochemically  active  electrode.  Only  a  few 
of  the  published  computational  studies  have  investigated  the 
effects  of  electrode  structures  on  cell  performance  [13-15].  These 
studies  are  based  on  volumetric  models  that  use  a  homogenous 
domain  characterized  by  macroscopic  properties  such  as  porosity 
and  specific  surface  area  to  represent  the  porous  electrode/elec¬ 
trolyte  medium.  Accordingly,  these  models  reveal  macroscopic 
trends  such  as  decreasing  porosity  resulting  in  decreased 
permeability,  bulk  diffusion  coefficients  and  increased  bulk 
conductivity. 

In  this  study,  the  pore-scale  model  is  used  to  analyze  the  effect 
of  the  three  electrode  structures  that  are  reconstructed  from  XCT 
images  as  shown  in  Fig.  2.  Flereafter,  it  is  implied  that  increasing 
porosity  yields  a  decrease  in  surface  area  and  vice  versa.  The  tested 
structures  are  simulated  with  flow  rates  that  correspond  to  Re 
numbers  of  order  1CT5  to  10  under  galvanostatic  discharge  at 
400  A  nrr2  with  a  specified  inlet  SOC  of  50%.  Shown  in  Fig.  8  are  the 
distributions  of  the  half-cell  voltage  and  corresponding  applied 
pressure  drop.  The  results  show  that  higher  voltages  are  delivered 
from  the  cell  for  denser  electrodes  at  each  flow  rate.  Additionally, 
the  marginal  trend  in  voltage  with  respect  to  flow  rate  observed  for 
the  93.2%  porosity  structure  from  Fig.  5  is  also  extended  to  the  less 
porous  morphologies.  Under  these  operating  conditions,  cell 
voltage  is  expected  to  saturate  above  a  flow  rate  of  Re  =  0.05  for  all 
three  structures. 

Fig.  9  shows  the  local  distributions  of  normalized  V5+  concen¬ 
tration  and  current  density  for  the  three  electrode  structures  where 
flow  rate  is  maintained  on  the  order  of  Re  =  0.05,  under  a  galva¬ 
nostatic  discharge  of  400  A  m-2  and  a  specified  inlet  SOC  of  50%. 
Note  that  for  this  study,  flow  rate  is  held  constant  such  that  the 
averaged  velocity  is  on  the  order  of  Re  =  0.05.  The  convection- 
dominated  concentration  field,  characterized  by  small  gradients, 
is  shown  in  Fig.  9(1).  From  Fig.  9(2)  it  can  be  seen  very  clearly  that 
decreased  porosity  will  facilitate  a  more  uniform  current  density 
distribution  as  a  result  of  an  increase  in  the  fiber  density  and 
effective  conductivity  of  the  solid  electrode.  To  clarify,  electrodes 
with  low  effective  conductivities  result  in  higher  reaction  currents 
near  the  current  collector,  which  occur  in  order  to  limit  the  amount 
of  electronic  current  in  the  electrode.  In  the  scatter  plot  histogram 
of  Fig.  9(3)  it  can  be  concluded  that  for  the  same  flow  rate, 
decreased  porosity  will  reduce  both  the  range  and  gradients  of 
concentration  and  current  density. 

The  corresponding  averaged  spatial  and  statistical  distribution 
of  concentration,  overpotential  and  current  density  for  the  three 
electrode  structures  are  shown  in  Fig.  10.  Greater  concentration 
gradients  are  observed  for  lower  porosity  structures,  as 
observed  from  the  decrease  in  the  mean,  variation,  and  range  in 
electrolytic  concentration  in  Fig.  10(a).  This  observation  can  be 
explained  by  the  fact  that  the  concentration  gradient  is  contin¬ 
gent  on  both  the  electrical  discharge  rate  and  the  mass  transfer 
conditions.  Since  the  same  current  density  is  applied  in  all  three 
cases,  the  reaction  rate  associated  with  the  discharge  rate  should 
remain  the  same  for  all  three  electrode  structures  in  accordance 
with  Faraday’s  law.  Increased  porosity  leads  to  increased  effective 
diffusion  coefficients  and  hydraulic  permeability  [13],  resulting  in 
better  mass  transfer  conditions  [15].  In  addition,  the  decrease  in 
surface  area  associated  with  porosity  increase  yields  reduced 
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Fig.  8.  Effect  of  porosity  and  flow  rate  on  (a)  half-cell  voltage  and  (b)  applied  pressure 
drop.  Simulations  are  performed  for  three  different  electrode  structures  subject  to 
galvanostatic  discharge  of  400  A  m-2  with  a  specified  inlet  SOC  of  50%. 


reaction,  thereby  decreasing  concentration  gradients.  Our 
observation  that  lower  porosity  increases  concentration  gradi¬ 
ents  is  in  contradiction  with  Shah  et  al.  [13]  and  in  agreement 
with  You  et  al.  [15]. 

As  the  electrode  structure  decreases  in  porosity,  more  surface 
area  is  available  for  electrochemical  reactions,  and  as  a  result 
lower  current  densities  are  required  to  yield  the  same  bulk  reac¬ 
tion  rate  as  required  by  Faraday’s  law.  Fig.  10(b)  and  Fig.  10(c) 
shows  that  the  resulting  overpotential  and  current  density  fields 
will  be  more  uniform  and  have  lower  overall  magnitudes.  From 
Fig.  10(b2)  and  Fig.  10(c2),  it  can  also  be  seen  that  the  local  polarity 
near  the  current  collector  is  greater  for  higher  porosity  electrode 
morphologies.  Since  the  effective  conductivity  of  the  solid  elec¬ 
trode  will  increase  with  increasing  fiber  density,  there  will  be  less 
resistance  in  the  solid  phase,  and  so  more  reactions  are  expected 
to  occur  farther  away  from  the  current  collector,  as  seen  in 
Fig.  9(1).  Additionally,  the  higher  local  current  densities  of  the 
more  porous  structures  is  expected  to  contribute  to  higher  local 
ohmic  losses,  which  overall  result  in  the  lower  half-cell  voltages  as 
observed  in  Fig.  8(a). 
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Fig.  9.  Distribution  of  local  transport  properties  shown  on  (1)  a  slice  plot  for  the  normalized  V5+  concentration  (C/Cin  vs+)  about  the  Y-midplane,  (2)  a  surface  plot  of  the  current 
density  (/)  on  the  3D  electrode  fibers,  and  (3)  a  scatter  plot  histogram  of  the  corresponding  J  and  C/CinV5+  at  every  applicable  control  volume  on  the  electrode  structure  (using 
32  x  32  bins).  Results  are  shown  for  different  electrode  structures  corresponding  to  porosities  of  (a)  0.932,  (b)  0.897,  (c)  0.845  subject  to  flow  rates  of  Re  =  0.05,  galvanostatic 
discharge  of  400  A  m-2  and  specified  inlet  SOC  of  50%. 


4.3.  Effect  of  SOC 

The  bulk  electrolytic  concentration  will  change  continuously 
during  the  operation  of  the  VRFB.  This  effect  is  studied  by 
prescribing  different  SOCs  at  the  inlet  of  the  flow  cell.  The  three 
electrode  structures  from  Fig.  2  are  subject  to  galvanostatic 
discharge  of  400  A  nrr2  and  a  constant  flow  rate  equivalent  to 
Re  =  0.5  for  a  range  of  SOCs  specified  at  the  inlet  between  10%  and 


90%,  so  as  to  avoid  the  effect  of  oxygen  and  hydrogen  evolution  at 
higher  polarizations  which  are  not  accounted  for  in  the  present 
model  [15]. 

From  Fig.  11(a)  it  is  evident  that  the  marginal  gain  in  half-cell 
voltage  is  the  same  for  all  SOCs,  and  that  higher  output  voltages 
are  acquired  at  higher  SOCs.  The  latter  observation  can  be 
explained  by  the  fact  that  higher  SOCs  lead  to  larger  Nernst 
potentials,  and  so  smaller  localized  liquid/solid  voltage  drops  are 
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Fig.  10.  Effects  of  porosity  on  averaged  distributions  of  (a)  normalized  V5+  concentration,  (b)  overpotential,  and  (c)  current  density  along  (1 )  the  Z-direction  and  (2)  the  X-direction. 
Statistical  information  for  each  parameter  in  the  entire  domain  is  shown  in  (3),  where  the  dots  indicate  the  mean  value,  the  capped  bars  indicate  the  variation  to  one  standard 
deviation,  and  the  uncapped  lines  indicate  the  minimum  and  maximum  values.  Results  are  shown  for  different  electrode  structures  corresponding  to  porosities  of  (a)  0.932,  (b) 
0.897,  (c)  0.845  subject  to  flow  rates  of  Re  =  0.05,  galvanostatic  discharge  of  400  A  m-2  and  specified  inlet  SOC  of  50%. 


sufficient  to  facilitate  electrochemical  reactions.  Due  to  these 
changes  in  Nernst  potentials,  it  is  important  to  note  that  inte¬ 
grated  average  overpotential  alone  is  not  sufficient  in  predicting 
output  voltages  [21].  For  instance,  the  symmetric  overpotential 
profile  about  SOC  =  50%  [8,15]  as  seen  in  Fig.  11(b)  may  mislead 
some  into  believing  that  both  high  and  low  SOCs  will  result  in 
higher  output  voltages. 

In  another  study,  a  single  electrode  structure  (the  93.2%  porosity 
morphology  in  Fig.  2(a))  is  subject  to  a  range  of  flow  rates  (Re  of 
order  10”5  to  10)  as  well  as  a  range  of  SOCs  (from  10%  to  90%)  under 
galvanostatic  discharge  of  400  A  m-2.  From  Fig.  12(a),  it  appears 


that  higher  SOCs  result  in  higher  overall  voltages  for  all  flow  rates. 
Note  that  in  the  extreme  case  of  low  flow  rate  (Re  <  0.0005)  and 
low  SOC  ( ~  10%),  a  dramatic  decrease  in  output  voltage  is  observed. 
This  observation  can  be  explained  by  inspecting  the  averaged 
Z-direction  concentration  profile  for  three  flow  rates  at  10%  SOC.  As 
shown  in  Fig.  12(b),  below  a  flow  rate  of  Re  =  0.0005,  a  fuel  star¬ 
vation  scenario  exists,  where  the  concentration  gradient  is  so  steep 
that  all  of  the  electrolytic  species  is  consumed  within  the  height  (Z) 
of  the  flow  cell.  It  can  be  seen  that  cell  performance  will  vary 
strongly  as  electrolytic  concentration  changes  during  operation  of 
the  VRFB. 
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Fig.  11.  Effects  of  porosity  on  (a)  half-cell  voltage  and  (b)  overpotential  as  functions  of 
specified  inlet  SOC  under  galvanostatic  discharge  of  400  A  m-2  and  constant  flow  rate 
of  Re  =  0.5. 


Fig.  12.  Effects  of  flow  rate  on  (a)  half-cell  voltage  for  different  specified  inlet  SOCs.  (b) 
Normalized  V5+  concentration  in  Z-direction  for  specified  inlet  SOC  of  10%.  Fuel  star¬ 
vation  is  observed  for  the  Re  =  0.0005  case.  Simulation  is  performed  for  the  93.2% 
porosity  structure  under  galvanostatic  discharge  400  A  m-2. 


5.  Conclusion 

An  isothermal  3D  pore-scale  model  [8]  is  used  to  investigate  the 
effects  of  electrolyte  flow  rate,  electrode  morphology,  and  inlet  SOC 
on  performance  characteristics  of  VRFBs.  Virtual  electrode  struc¬ 
tures  are  reconstructed  from  XCT  data  of  a  porous  carbon  felt 
material  and  analyzed  in  3D  using  the  pore-scale  model.  Localized, 
averaged  and  statistical  distributions  of  concentration,  over¬ 
potential  and  current  density  fields  across  the  tested  3D  electrode 
domains  are  determined  to  identify  the  key  parameters  affecting 
the  electrolyte  flow  and  utilization.  The  results  show  that  overall 
cell  voltage  increases  with  increasing  flow  rate  for  the  investigated 
electrode  morphologies  and  SOCs  as  a  result  of  decreasing 
concentration  gradients.  However,  the  marginal  gain  in  cell  voltage 
diminishes  once  the  concentration  field  approaches  uniformity  at 
sufficiently  high  flow  rates.  It  is  also  observed  that  cell  performance 
improves  for  denser  electrode  structures  with  more  active  surface 
area,  as  a  result  of  more  uniform  and  lower  absolute  current  density 
and  overpotential  fields.  Finally,  the  simulations  suggest  that 
a  significant  detriment  to  the  performance  of  a  VRFB  occurs  in  the 


event  of  fuel  starvation  for  low  electrolytic  concentrations  and  low 
flow  rates. 
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Nomenclature 

C:  concentration  [mol  m-3] 

D:  diffusivity  [m2  s-1] 

E°:  open  circuit  voltage  [V] 

F:  Faraday’s  constant  [C  mol-1] 

H:  total  cell  height  [m] 

J:  surface  current  density  [A  m-2] 


k:  reaction  rate  constant  [m  s-1] 

I:  average  pore  size  [pm],  total  cell  length  [m] 
h:  surface  unit  normal 
N:  mole  flux  [mol  m-2] 

P:  fluid  pressure  [Pa] 

R:  universal  gas  constant  [J  mol-1  K-1] 

Re:  Reynold’s  number 

SOC:  state  of  charge 

T:  operating  temperature  [K] 

u:  velocity  of  the  electrolyte  flow  [m  s-1] 

W:  total  cell  width  [m] 

X:  component  in  the  X  direction  [mm] 

Y:  component  in  the  Y  direction  [mm] 

Z:  component  in  the  Z  direction  [mm] 

Greek 

a:  transfer  coefficient 
p:  overpotential  [V] 

0:  potential  [V] 

k:  conductivity  [S  m-1] 

v:  electrolyte  ldnematic  viscosity  [m2  s-1] 

p:  electrolyte  density  [kg  m-3] 

Subscript 

a:  anodic  reaction  quantity 
avg:  average  quantity 
c:  cathodic  reaction  quantity 

ext:  externally  applied  quantity,  current  collector  quantity 
e:  liquid  electrolyte  property 
j:  species  CjJe  {V(IV),  V(V),  H+,  H20,  SO^“} 
s:  solid  fiber  phase  property 

Superscript 

*:  normalized  quantity 
0:  equilibrium  or  initial  quantity 
in:  cell  inlet  value 
s:  surface  property 
+:  positive  half-cell  quantity 
— ;  negative  half-cell  quantity 


